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ABSTRACT 

Aims. Due to the importance of accurate Fourier parameters, we devise a method that is more appropriate for deriving these parameters on 
low-quality data than the traditional Fourier fitting. 

Methods. Based on the accurate light curves of 248 fundamental mode RR Lyrae stars, we test the power of a full-fetched implementation of 
the template method in the computation of the Fourier decomposition. The applicability of the method is demonstrated also on datasets of filter 
passbands different from that of the template set. 

Results. We examine in more detail the question of the estimation of Fourier-based iron abundance [Fe/H] and average brightness. We get, for 
example, for light curves sampled randomly in 30 data points with a = 0.03 mag observational noise that optimized direct Fourier fits yield 
cr([Fe/H]) = 0.33, whereas the template fits result in cr([Fe/H]) = 0.18. Tests made on the RR Lyrae database of the Large Magellanic Cloud 
(LMC) of the Optical Gravitational Lensing Experiment (OGLE) support the applicability of the method on real photometric time series. These 
tests also show that the dominant part of error in estimating the average brightness comes from other sources, most probably from crowding 
effects, even for under-sampled light curves. 
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1. Introduction 

The method of Template Fitting (TF) is a widely used approach 
in data analysis. In astronomical applications we find examples 
from spectrum analysis (e.g., Bertone et al. 2004) to galaxy 
classification and redshift estimation (e.g., Wolf, Meisenheimer 
& Roser 2001; Padmanabhan et al. 2005). The method is based 
on the simple assumption that the population to be studied 
contains targets sharing the same topological properties as the 
members of the template set. The latter is defined as a set 
containing all possible 'flavors' of the given population and 
possessing accurately known parameter arrays - e.g., spec- 
trum, redshift, etc. The actual implementation of the TF method 
ranges from the simple few-template direct match (e.g., Jones, 
Carney & Fulbright 1996; Layden 1998) to the more sophis- 
ticated Artificial Neural Network method (e.g., Collister & 
Lahav 2004). 

Here we present a 'brute force' direct TF method that is 
aimed at the computation of the Fourier decompositions of fun- 
damental mode RR Lyrae (RRab) stars. The goal of this inves- 
tigation is to provide a reliable method for the computation of 



the Fourier decomposition of any observed RRab light curve 
even if observational noise or poor sampling impair standard 
Fourier decomposition. Our approach is different from that of 
Kanbur & Mariani (2004) and Tanvir et al. (2005), who em- 
ployed principal component analysis (PCA) to parametrize the 
light curves of RR Lyrae and Cepheid variables. With the aid of 
PCA one is able to create a smaller set from the templates and 
represent the target by a low-degree (in some sense optimum) 
PCA decomposition. However, low-order PCA decompositions 
are often insufficient if more subtle features are required to fit 
(see also Fig. 2 of Tanvir et al. 2005). In addition, in the case 
of RR Lyrae stars, Blazhko effect further increases the possible 
types of light curves(see also Jurcsik, Benko & Szeidl, 2002). 
Therefore, we are resorted to a method that is able to handle a 
large variety of light curves, flexible enough but does not 'over- 
fit' the data. 

Compared to earlier related works on RR Lyrae stars, the 
present one utilizes a much larger template set, containing 248 
accurate RRab light curves. The method is tested through a 
comparison with an optimized Fourier fit. We focus on the ac- 
curacy of the estimation of the iron abundance [Fe/H] based 
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on the Fourier decomposition (Jurcsik & Kovacs 1996, here- 
after JK96, see also Kovacs 2005) and on the determination of 
the period-luminosity-color (PLC) relation (Kovacs & Walker 
2001, hereafter KW01). Neither of these quantities can be 
accurately estimated with direct Fourier fit if the number of 
data points is low and the noise is high, such as in the case 
of the V and B band observations of OGLE (Soszynski et 
al. 2003). Furthermore, the success of any utilization of the 
color indices in the computation of the physical parameters 
(most importantly that of T e s) strongly depends on the ac- 
curacy by which we estimate average colors. As far as the 
computation of the Fourier-based [Fe/H] is concerned, here 
the accurate estimation of <pn is also important, because of 
the strong dependence of [Fe/H] on this quantity - see the 
application of the empirical [Fe/H] formula on the MACHO 
LMC data by Kunder et al. (2006). Current efforts in de- 
riving [Fe/H] on large samples of stars in globular clusters 
and nearby galaxies from low/medium-dispersion spectroscopy 
(e.g., Sandstrom, Pilachowski & Saha 2001 [M3]; Gratton et al. 
2004 [LMC]; Clementini et al. 2005a [NGC 6441]; Clementini 
et al. 2005b [Sculptor dwarf spheroidal galaxy]; Sollima et al. 
2006 [cj Cen]) make the accurate computation of the Fourier 
decompositions of RRab stars even more interesting. 

In the subsequent sections we describe the method, opti- 
mize the template function, discuss template completeness, in- 
vestigate the effect of choosing different filter passbands for 
the target and template sets and test the accuracy of [Fe/H] and 
average magnitudes derived from the TF method. Finally, we 
present results based on a limited dataset from the OGLE RRab 
database. 



- The unbiased estimate of the fitting accuracy (the r.m.s. of 
the residuals between the fit and the data) is minimum. 

- The total amplitude of the fitted curve is not greater than 
fjA m , where A, ot is the total amplitude of the target sig- 
nal and fj is our second empirical parameter, and is also 
chosen to be equal to 1 .2. 
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Fig. 1. Amplitude ratios of the successive Fourier components 
of the 248 RRab stars of the basic dataset. For better visibility, 
data points related to the individual objects are shifted in the 
horizontal direction. Except for orders above 9-10, for most 
of the objects the amplitudes monotonically decrease with the 
increase of the order of the components. 



2. The TF and the direct Fourier methods 

In order to make a meaningful (and fair) comparison between 
the Fourier parameters derived from the TF method and those 
obtained from a straightforward Fourier fit, we need to perform 
the latter in an 'optimum way' (i.e., the way it would be done by 
a skilled data-analyst by checking the fit for different orders and 
avoiding - if possible - overshooting and appearance of strong 
wiggles for under-sampled light curves). First we describe our 
automated method for the Direct Fourier Fitting (DFF) and then 
give details on the Template Fourier Fitting (TFF 1 ). 

In the case of DFF we tried to establish a set of criteria 
that results in 'good looking' fits and can be applied automat- 
ically, without further human inspection. Because we employ 
standard unweighted least squares Fourier fits, the only param- 
eter by which we can influence the quality of the fit is the order 
of the Fourier sum. By scanning the orders from 1 to 10, we 
choose the highest order at which the following criteria are sat- 
isfied: 

- The Fourier amplitudes are still 'nearly' monotonically de- 
creasing (i.e. A;_i > /aAi, where f A is one of our 'trial and 
error' parameters and is set equal to 1 .2, based on the in- 
spection of the observed Fourier amplitudes of RRab stars 
- see Fig. 1). 

1 A Fortran'77 source code of TFF is available at 
http://www.konkoly.hu/staff/kovacs/tff.html 



These criteria make DFF reasonably stable for poorly sam- 
pled light curves, and produce accurate fit for well-sampled 
ones. We note that Ngeow et al. (2003) employed a somewhat 
similar method in deriving smooth and stable Fourier decom- 
positions for Cepheid light curves. Their method employs 'sim- 
ulated annealing' (see Press et al. 1992) constrained by the pe- 
riod dependence of the Fourier amplitudes of Cepheids (the so- 
called Hertzsprung progression). 

For TFF, our approach is similar to that of Layden (1998), 
except that: (i) we use a much larger template set, based on 
individual variables and not on a limited set of visually selected 
classes; (ii) we allow low-degree polynomial transformation of 
the template in finding the best fit. 

First we choose a set of Fourier decompositions derived 
from well-observed, densely sampled light curves. Such a set 
is available from our earlier studies on the V band light curves 
of RRab stars (see the CDS archive of KW01). There are alto- 
gether 492 variables, with 105 stars from the Galactic field and 
the rest from various globular clusters and from the Sculptor 
dwarf galaxy. Because some of variables are poorly sampled, 
we apply the following selection criterion in order to employ 
only the best quality light curves. First we define the Quality 
Factor (QF) of a given light curve as a quantity proportional to 
the ratio of the amplitude A \ of the first Fourier component to 
its simple error estimate, i.e. 

QF = y/NAx/o-m , (1) 
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where N is the number of data points and cr fit is the unbiased 
estimate of the standard deviation of the residuals of the fit. The 
quality control is established by requiring QF > QF min , where 
QF min is a preset threshold. By changing QF m i n from 100 to 
200, the number of the remaining stars decreases from 336 to 
193. Finally we decided to apply QF m ; n = 150 and obtained 
a set of 251 stars. The overwhelming majority of these stars 
have N > 100, and there are only three stars with N < 40. 
Since at lower number of data points there is a greater risk of 
having erroneous Fourier decompositions, we omit these three 
stars (FH Vul, IV Hya and NGC1841 V4). Finally we arrive to 
our basic dataset containing 248 variables. We use this set as 
the template set throughout this paper. 

Once the template set is selected, for a target light curve we 
find the best fitting template in the following way. 

- Compute densely sampled folded light curves from the 
template Fourier decompositions. We denote these func- 
tions by {xij(if)}, where subscript i stands for the array in- 
dex in the folded light curve and j refers to the template 
identification. Phase <p of the template is arbitrary at this 
step. 

- Compute folded light curve {T,} for the target. 

- For each initial phase <p and for each template, minimize 
the following quantity: 



1 N 

1=1 



(2) 



where 



C + Xij((f) 



AM 



if M = 
if M> 1 



Here M denotes a preset polynomial degree, to be deter- 
mined in Sect. 3 as a data quality-dependent parameter. 
While scanning <p, we employ quadratic interpolation for 
the template in order to get a good approximation for its 
value at the moments where the target is given. 

- TFF is computed by the Fourier decomposition of that 
{Xi(ip)}, which minimizes 2)j(tp) of Eq. (2). 

- The above steps are to be supplemented by the ones to be 
discussed in Sect. 3 for the optimum choice of M. 

Typically we use template light curves sampled in 300 points 
and require an accuracy of 10~ 5 in the phase match between 
the target and template. The search for optimum phase is made 
iteratively, starting with some 50 phase steps. Execution time 
is not an issue with current several GHz machines. 

One can construct other types of TF methods, by using 
different functional dependence of the target on the template 
members (e.g., linear or polynomial multi-template functions). 
However, in our approach we consider the current template set 
only as a subset of an ever growing master set that will be ac- 
cumulated in the future. For this ideal set we might need only 
a scaling factor for a very precise fit of any target, because the 
master set will contain all 'flavors' of RRab stars and the fit- 
ting routine needs to perform the search only among the single 
template members. In addition, more complicated functional 
dependence would make the fitting procedure slower. 



3. Template polynomial degree and completeness 

Before testing the TFF method described in Sect. 2, it is neces- 
sary to determine the best polynomial degree M to be used in 
the template transformation. Furthermore, it is also important 
to examine if the template set with the adopted TFF method is 
capable of reproducing all (or most of) the light curve 'flavors' 
observed among RRab stars. This latter property is connected 
to what we call 'completeness' and to be defined somewhat 
more precisely later in this section. 

The purpose of introducing the polynomial template trans- 
formation is to increase our freedom in reaching higher accu- 
racy in fitting targets. Obviously, employing a too high-degree 
polynomial may lead to instability, similarly as if we used high- 
order Fourier fit. Because the prime goal of the application of 
the template method is just to avoid this type of instability, we 
accept the lowest polynomial degree that yields fits of similar 
quality as the higher degree ones. 

In order to rank the results obtained with various polyno- 
mial degrees, we need to define a function that characterizes the 
quality of the fit for an ensemble of targets. Using some aver- 
age of the standard deviations of the fits to the individual targets 
is not satisfactory, because then, poorly fitted small-amplitude 
variables may stay hidden due to the small standard deviations 
associated with their small amplitudes. A possible normaliza- 
tion by the amplitude (see Eq. (1) for QF) may partially cure 
this problem, but we found more satisfactory to use a function 
that is independent of the amplitude and more closely related 
to our prime interest in deriving accurate phases. We introduce 
the following quantity to characterize the goodness of the fit 



RMS(A<£) 



4 



3 2j^ A ^i 



k=2 



(3) 



where Aifki is the difference between the target and the best 
template-fitted phases. The epoch-independent phase is defined 
in the usual way <pk\ = <fik - k<p\ (Simon & Lee 1981). Notice 
that the above expression utilizes all three low-order phases, 
not only ip^ that enters in the empirical formula for [Fe/H]. 
We think that using more phases makes the results more sta- 
ble against statistical fluctuations. At the same time, adding 
high-order phases would make the above statistic biased toward 
small fitting errors, that we would like to avoid. 

By using RMS(A^), the optimum polynomial degree M is 
determined in the following way. For any fixed M we start with 
the 248 stars of the basic dataset and for each variable of this 
set we find the best fitting template selected from the remain- 
ing 247 stars. The time series of the target is computed from the 
Fourier series given in the basic dataset for the target chosen. 
The time base of the sampling is taken arbitrarily as 12.3456P, 
where P is the period of the target. The sampling rate is quasi- 
uniform with a small randomness of the size of the exact uni- 
form sampling. Tests are run both with and without noise added 
to the synthetic data. We note that other choices of the time 
base would also serve the purpose, except for near integer mul- 
tiple of the period, when the chance of regular sampling of the 
phased light curve would be greater. 

Once the best template (i.e., {Xi(<p)}, in Eq. (2)) is found, 
the resulting Fourier decomposition is compared with that of 
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the target via Eq. (3). In this way at each fixed M we get 248 
RMS(Ai/j) values that can be analyzed statistically and compare 
with those obtained at other M values. 

The most straightforward way to compare the RMS(Ai/j) 
values is to compute their probability distribution function 
(PDF). In Fig. 2 we show such functions computed for the 
noiseless (i.e., <x = 0.0) simulations with N = 20. As expected, 
the distribution functions become more narrow and concen- 
trated at lower RMS(A^) values when M is chosen to be low. 
From this test the optimum M at N = 20 is expected between 2 
and 4 for noiseless data. 
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Fig. 2. Probability distribution functions of the average phase 
differences RMS(A</j) (Eq. 3) for test targets at various tem- 
plate polynomial degrees M. All target light curves were gen- 
erated with N = 20 data points. No noise is added. The test was 
performed on the basic dataset. Each PDF is normalized to its 
maximum. We draw attention to the narrowing of the PDFs at 
low M values of 2^4. 

If we perform the above test at other values of N, then we 
get different values for the optimum M. This is understandable, 
because at low the high polynomial degree leads to stronger 
instabilities, whereas at high we are able to fit templates of 
high degree, thereby reaching higher accuracy. In order to get 
an estimate on the size of this shift of M, we performed addi- 
tional tests with N = 15, N = 40 and N = 100. For an easier 
comparison, instead of plotting PDFs, we compute the num- 
ber of stars that have RMS(A</)) greater than RMS max , where 
the latter quantity is chosen in a way which ensures that, at 
least for some distribution functions, the number of stars sat- 



isfying this condition is small. From the inspection of Fig. 2 
we choose RMS max = 0.12, because for weakly-spread distri- 
bution functions the tail seems to be separated from the bulk 
of the distribution at this value. (We note that our conclusion 
does not change by choosing other cutoff values in the range of 
0.07-0.17.) 
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M 

Fig. 3. Dependence of the relative number of poorly fitted tar- 
gets on the polynomial template degree M and on the number 
of data points Af of the target time series. Tests were made on 
the basic dataset with N tot = 248 stars. No noise was added to 
the synthetic data. 
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Fig. 4. As in Fig. 3, but with added noise as given in the shaded 
box. 

Fig. 3 shows that the size of the shift in the location of the 
minima of the functions is small for low Af values and we can 
still stay in M — 2-4, with hitting the low and high boundaries 
at low and high Af values, respectively. For modest (or high) Af 
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values the minima become shallower, therefore, there will not 
be much difference between choosing moderately low M and 
the optimum one. However, it is clear that at large N the mini- 
mum will be lower and gradually shifted to rather high M val- 
ues. On the other hand, these high-A^, well-sampled cases are 
not interesting from the present point of view, because these are 
well- treated by standard Fourier fitting methods. (Exceptions 
are of course cases when high noise prohibits the traditional 
approach, and we will become better off again by using the 
template method at low M.) 

To check the effect of noise on the optimum polynomial de- 
gree, we repeated the above test by adding moderate Gaussian 
noise of <x = 0.03 mag to the synthetic target signals. The re- 
sult is shown in Fig. 4. As expected, the properties observed at 
low N in the noiseless case of Fig. 3 are shifted to higher N. At 
the same time, the low-Af signals get basically out of control at 
high M. 

Due to the dependence of the optimum template polyno- 
mial degree on data quality, we need to examine this relation 
more closely. The data quality is characterized by the Signal- 
to-Noise Ratio (SNR), defined similarly to QF 



SNR = VNA/cr fit , 



(4) 



where A is the total amplitude of the light variation, cr fit is the 
unbiased estimate of the standard deviation of the residuals be- 
tween the target and the fit. We performed tests similar to the 
ones described above and examined the behavior of RMS(A^) 
as a function of SNR and M. The SNR values were obtained 
from the results at M — 1 . For N, M and cr we took the fol- 
lowing values: N = 15, 30, 45, 60, 90, M = 0, 1, 2, 3, 4 
and cr = 0.0, 0.03, 0.06, 0.09. Simulations corresponding to 
the same M were averaged and plotted as functions of SNR 
in Fig. 5. The errors of the averages are about the size of the 
circles, except toward the low SNR end, where they increase 
to 0. 1-0.2. It is clear from this figure that the optimum polyno- 
mial degree is a function of SNR. Therefore, the TFF algorithm 
described in Sect. 2 is supplemented by the following steps in 
selecting the optimum polynomial degree. 

- Compute SNR by applying TFF at M = 1 

- Choose the best M depending on the computed SNR 

f 0, if SNR < 50.0 
M = \ 1, if 50.0 ^ SNR ^ 150.0 
( 2, if SNR > 150.0 

- Compute TFF with M determined above 

We note that SNR can also be estimated by the simplest M = 
fit, but we choose M = 1 , because it yields a somewhat more 
consistent result (i.e., cleaner separation of the averages corre- 
sponding to the various M values) and because of the unavoid- 
able scaling in fitting light curves of different passbands (see 
Sect. 5). In the rest of this paper we use TFF with the above 
optimized M. 

Next we address the question of template completeness. As 
already mentioned, we would like to measure the ability of the 
basic dataset in reproducing each member of the set with the 
aid of the TFF method. It is clear that we need to define the 
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Fig. 5. Dependence of the average of RMS(A^) on SNR at var- 
ious template polynomial degree M shown in the inset. The 
figure is based on simulated data described in the text in detail. 

meaning of the word 'reproduce'. Here, employing the same 
argument as earlier, we use RMS(A^) as a quantity characteriz- 
ing the goodness of the fit. Then we say that the template set is 
complete at the level of Ctff at RMS max , if the fraction of vari- 
ables that satisfy the condition RMS(A</>) < RMS max is equal 
to Ctff- Because the introduction of Ctff is aimed at to char- 
acterize the ability of the template set to 'reproduce' itself by 
using TFF, we generate noiseless, well-sampled, high-Af light 
curves to derive the dataset necessary for the computation of 
Cjff- The result of is shown in Fig. 6. We see that the com- 
pleteness is close to 90% at RMS max = 0.1 and reaches 95% 
at RMS max = 0.15. For some 80% of the stars we get matches 
with smaller RMS than 0.05. 




RMS„ 



Fig. 6. Completeness function of the TFF method. The sym- 
bol Ctff denotes the relative number of variables that have 
RMS(A^)) < RMS max . We used the 248 light curves of the ba- 
sic dataset with N = 200 synthetic data points per object and 
polynomial degree M — 2. 



Although the above numbers indicate a reasonable com- 
pleteness, we note that there are stars that stubbornly resist to 
accurate template fitting and prevent high completeness even at 
low-accuracy (e.g. for RMS max > 0.15). For example, variable 
M107 V12, cannot be fitted, yielding RMS(Ay) = 0.263 and 
A^3i = i/?3i (target) - tpn (template) = 0.389. There are some 
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15 stars that have \A<pn\ > 0.1. The fact that some stars cannot 
be fitted even if their light curves are very densely sampled, 
can be explained by one (or some) of the following reasons: 
(i) the size of the current template set is small, and therefore, 
it is unable to reproduce some of the existing light curves with 
a desirable accuracy; (ii) some of the stars may exhibit long- 
periodic amplitude- and phase-modulations that result in dis- 
crepant template matches; (iii) the data on the targets were in- 
sufficient (too few data points, gaps in the folded light curves, 
etc.), that has led to inaccurate Fourier decompositions (in spite 
of these stars passing through our criteria for template mem- 
bership); (iv) instrumental effects (e.g. daily or seasonal drifts) 
make the given light curve unique and therefore, not treatable 
by the template method; (v) some stars have periods close to 
integer ratios of one day, that again, may lead to unique light 
curves, especially if it is combined with property (iv). We think 
that from the present dataset we cannot decide which of the 
above effects is responsible for the outlier status of some of the 
stars (a closer examination of some of the outliers has shown 
that one can find examples/suspects for all these five possibili- 
ties). Considering that the current template set contains several 
stars with limited coverage either in time or in amount of data 
points, the existence of the few outliers is not surprising. It is 
clear that the present template set should be further extended 
by utilizing new observations and by leaving out objects with 
poor quality light curves. 

To illustrate the difference between the matches produced 
by DFF and TFF at low number of data points, in Fig. 7 we 
show two examples. Although for other data distributions DFF 
might behave less erratically, the behavior shown is quite com- 
mon at low number of data points. We note that in many cases 
the best matching template may have widely different period 
from that of the target. For instance, in the example shown, DX 
Del has a period of 0.47262 d, whereas the best match is pro- 
duced by u Cen V176, that has a period of 0.74275 d. For RR 
Leo the situation is different, here we have P = 0.45239 d for 
the target and P = 0.45930 d for the matching template SW 
Aqr. 

4. Estimation of the photometric [Fe/H] 

From the point of view of applications, it is important to ex- 
amine the accuracy of the TFF method in determining the 
photometric [Fe/H] with the aid of the formula of JK96. We 
used the basic dataset in generating simulated light curves with 
the following parameter values: N = 15, 30, 45, 60, 90; 
cr = 0.0, 0.03, 0.06. The polynomial degree M was opti- 
mized as described in Sect. 3. For each N, cr and for each star, 
we compute A[Fe/H], the difference between the [Fe/H] val- 
ues computed from the target and from the fitted light curves. 
From these values we derive: <x(A[Fe/H]), the standard devia- 
tion of the [Fe/H] differences; R01, the ratio of the stars with 
|A[Fe/H]| < 0.1 to the total number of stars (i.e. to 248); T/D, 
the ratio of the number of stars for which |A[Fe/H]| is smaller 
for TFF than for DFF to the ones for which the opposite is true. 

In Table 1 we show the average values obtained for these 
quantities for the two methods. The following conclusions can 
be drawn from this table: 
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Fig. 7. Examples of the performance of TFF at low data point 
numbers. Dots, dashed and continuous lines are for the noise- 
less synthetic data, DFF and TFF fits, respectively. The noise- 
less synthetic data have been obtained by the sparse sampling 
of the accurate observed light curve. 

- Except at low noise level and high (i.e. £ 45) data point 
number, TFF always yields more accurate [Fe/H] values in 
the average sense. 

- Within the above limit, the number of accurate [Fe/H] esti- 
mates (i.e. those with |A[Fe/H]| < 0.1) is always larger for 
TFF. 

- Within the above limit, the number of cases when TFF 
yields smaller error than DFF is always larger than that of 
the opposite situation. 

- This better performance is especially well visible at higher 
noise levels. 
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Fig. 8. Standard deviation of the difference between the target 
and computed Fourier [Fe/H] values as a function of SNR. The 
test was made on the 248 stars of the basic dataset. Error bars 
show the +lcr ranges of the means. 
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Table 1. Accuracy of the determination of [Fe/H]. 
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4.4 


0.00 


30 


0.21 


0.11 


0.7 


0.9 


1.5 


0.00 


45 


0.12 


0.09 


0.9 


0.9 


0.6 


0.00 


60 


0.03 


0.09 


1.0 


0.9 


0.3 


0.00 


90 


0.02 


0.08 


1.0 


0.9 


0.2 


0.03 


15 


0.76 


0.27 


0.2 


0.4 


3.9 


0.03 


30 


0.33 


0.18 


0.4 


0.5 


1.9 


0.03 


45 


0.24 


0.17 


0.5 


0.6 


1.4 


0.03 


60 


0.13 


0.15 


0.7 


0.7 


1.0 


0.03 


90 


0.12 


0.13 


0.7 


0.7 


0.8 


0.06 


15 


0.86 


0.37 


0.1 


0.3 


4.9 


0.06 


30 


0.45 


0.29 


0.3 


0.3 


1.6 


0.06 


45 


0.36 


0.22 


0.3 


0.4 


1.9 


0.06 


60 


0.24 


0.19 


0.4 


0.5 


1.8 


0.06 


90 


0.23 


0.18 


0.4 


0.5 


1.5 



Note: <x: standard deviation of the noise added to the synthetic light 
curves that were generated from the template Fourier decompositions; 
N: number of data points; cr D : standard deviation of A[Fe/H] D = 
[Fe/H] t t - [Fe/H] DFF ; cr T : as cr D , but for TFF; R0l D : number of 
stars with |A[Fe/H] D | < 0.1 divided by the total number of stars 
of 248; R0l T : as R0l D but for TFF; T/D: number of stars with 
|A[Fe/H] T | < |A[Fe/H] D | divided by the number of stars satisfying the 
opposite inequality. The result is based on the 248 stars of the basic 
dataset. 

To make a comparison yet in another parameter domain, in 
Fig. 8 we plot the standard deviation of A[Fe/H] as a function 
of SNR. Because the various simulations yield different indi- 
vidual SNR values, the total range of SNR was divided into 
20 bins and the standard deviations of the various A[Fe/H] val- 
ues within these bins have been computed. Except for low SNR 
values (i.e. for log SNR < 1.75), all bins contain some 100-300 
simulations (at log SNR — 1.5 we have only 10). Although the 
scatter within the bins is very large, the averages are fairly accu- 
rately estimated and the difference between the two methods is 
clearly visible. This figure may give some guidance to a rough 
error estimation of the methods. In general, we may expect 
rather large errors - in the average sense - if log SNR < 2.0. 
The average errors can be substantially decreased for DFF for 
log SNR > 2.0 if we employ 3cr clipping on the derived [Fe/H] 
values. In this way the two methods will perform nearly in the 
same way for log SNR > 2.0. The effect of 3cr clipping on TFF 
is minimal, because the number of outliers is much lower when 
TFF is employed. As an example, at log SNR = 2.75, with 
3<x clipping we loose 23% of the stars in the case of the DFF 
method. The same figure for TFF is only 2%. 

5. Estimation of the average magnitudes in 
various colors 

We test the applicability of the TFF method in estimating the 
average magnitudes of light curves observed in various wave- 
bands. Unfortunately, the amount of good-quality multicolor 
data available for us is much lower than that of the single-color 
data. For simplicity, we take a set from the globular cluster data 
used by KW01 for the derivation of the PLC relation for RRab 



stars. The present set contains B, V and I light curves from 
globular clusters NGC 1851, 4499, 6362 and 6981. There are 
95, 95 and 83 variables in B, V and in I colors, respectively. 
Although the amount and the quality of these data are less fa- 
vorable than the ones used for testing the accuracy of the pho- 
tometric [Fe/H], they are sufficient to get rough error estimates 
on the determination of average magnitudes by using different 
methods. 
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Fig. 9. Empirical probability distribution functions of 
RMS(A</j) for noiseless test signals drawn from samples 
of RRab stars of different wavebands. All test signals have 
N — 90 data points. The TFF method with the basic dataset 
(available only in color V) is used in the computation of Atp. 

Before discussing results on the average magnitudes, it is 
interesting to compare the accuracy of the fit in the different 
colors. As before, we use RMS(A^) to characterize the good- 
ness of fit. Figure 9 shows the PDFs of this quantity for the 
different colors. As expected, the V-band data are fitted most 
accurately. The B and I data perform similarly, with a slight 
preference (but probably within the error limits of the present 
test) toward the B data. The fact that these, seemingly very sim- 
ilar light curves are not transformable by TFF, shows that they 
contain independent pieces of information on the pulsation. In 
the present context this result suggests to avoid the use of tem- 
plates of different waveband from that of the target if we are 
aimed at the estimation of the Fourier decomposition with the 
aid of the TFF method. 

In comparing the average magnitudes, we proceed in the 
same way as in the test of [Fe/H] in Sect. 4, except that now the 
target sets are limited on the cluster data mentioned above. The 
results are ranked on the basis of the standard deviations com- 
puted from the differences between the averages of the target 
(as given by the zero frequency constant in its accurate Fourier 
decomposition) and the estimated values. Table 2 shows the re- 
sult of the computation for the three colors. 

Although other realizations (determining the data distribu- 
tion and noise) lead to somewhat different results, the following 
basic trends seen in the table survive. 

- The simple arithmetic average (AVE) has the largest scatter 
and it is almost independent of the noise level. 
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Table 2. Accuracy of the magnitude averages. 



Band 


AV/N 15 


30 


45 


60 


90 


cr = 0.0 


B: 


AVE 


0.082 


0.051 


0.032 


0.027 


0.017 




DFF 


0.037 


0.009 


0.004 


0.002 


0.001 




TFF 


0.030 


0.007 


0.005 


0.003 


0.002 


V: 


AVE 


0.063 


0.040 


0.025 


0.021 


0.013 




DFF 


0.032 


0.007 


0.004 


0.002 


0.001 




TFF 


0.010 


0.005 


0.003 


0.002 


0.002 


I: 


AVE 


0.039 


0.023 


0.013 


0.013 


0.008 




DFF 


0.023 


0.004 


0.003 


0.002 


0.001 




TFF 


0.012 


0.004 


0.002 


0.001 


0.001 


cr = 0.03 


B: 


AVE 


0.082 


0.051 


0.032 


0.027 


0.018 




DFF 


0.043 


0.011 


0.007 


0.005 


0.004 




TFF 


0.027 


0.008 


0.007 


0.005 


0.004 


V: 


AVE 


0.064 


0.040 


0.024 


0.022 


0.014 




DFF 


0.030 


0.014 


0.006 


0.005 


0.003 




TFF 


0.022 


0.007 


0.006 


0.004 


0.003 


I: 


AVE 


0.040 


0.024 


0.013 


0.014 


0.008 




DFF 


0.024 


0.008 


0.006 


0.005 


0.004 




TFF 


0.015 


0.007 


0.005 


0.004 


0.003 


cr = 0.06 


B: 


AVE 


0.083 


0.052 


0.032 


0.028 


0.019 




DFF 


0.045 


0.020 


0.011 


0.009 


0.007 




TFF 


0.032 


0.013 


0.011 


0.008 


0.007 


V: 


AVE 


0.065 


0.041 


0.025 


0.023 


0.015 




DFF 


0.037 


0.018 


0.011 


0.008 


0.007 




TFF 


0.024 


0.012 


0.010 


0.008 


0.006 


I: 


AVE 


0.042 


0.027 


0.014 


0.015 


0.010 




DFF 


0.029 


0.015 


0.010 


0.008 


0.007 




TFF 


0.029 


0.013 


0.009 


0.008 


0.006 



Note: Each item of data corresponds to the standard deviation of the 
AA values, as given by the difference between the average magnitude 
of the target and the one obtained by the application of the various 
methods (AVE denotes the simple arithmetic average). The tests are 
based on cluster RRab stars as given in the text. The result depends 
somewhat on the realization used, but this does not change the basic 
trends shown here. 

- Within the statistical limits, TFF always performs better 
than DFF. 

- If the number of data points is greater than ~ 40, TFF and 
DFF yield results of similar accuracy. 

The near independence of AVE on the noise level is due to 
the fact that in the course of simple averaging the main source 
of error is the uneven distribution of data points in the phased 
light curve. It is seen that this effect is not negligible even at 
high number of data points. Therefore, for the accurate compu- 
tation of the averages, simple averaging should never be used. 
Although Fourier decompositions are less accurately estimated 
for light curves of colors different from that of the template set, 
averages have nearly the same accuracy in all colors. 

6. PLC estimates on artificial data 

The correlation of the period and the reddening-free magnitude 
Wb-v = V - Ry(B - V) is an important relation in estimat- 



ing RR Lyrae distances from colors in the visual wavebands 
(Dickens & Saunders 1965; Kovacs & Jurcsik 1997; KW01; 
see also Kovacs 2003 for application in Baade-Wesselink anal- 
ysis and Di Criscienzo, Marconi & Caputo 2004 for theoretical 
interpretation). Because of the large value of the selective ab- 
sorption coefficient Ry, random errors in the estimated mean 
magnitudes become amplified in Wb-v and thereby impair the 
accuracy by which we can employ the relation, among others, 
for distance determination. Therefore, it is worthwhile to test 
the effect of TFF in decreasing the error of Wb-v- 




-.5 -.4 -.3 -.2 -.1 -.5 -.4 -.3 -.2 



logio P logio P 

Fig. 10. PLC relations computed from tests on the 248 stars 
of the basic dataset. Test light curves were generated with 30 
data points and additive Gaussian noise of 0.03 mag standard 
deviation. The selective absorption coefficient Ry is set equal 
to 3.1. 

The test utilizes the 248 stars of the basic dataset. For each 
star we generate test light curves in the following way: 

- By using the Fourier decompositions of the V light curves, 
compute magnitude-averaged V and B colors from Eqs. (5) 
and (6) ofKWOl. 

- Compute zero-averaged synthetic light curves from the 
Fourier decompositions and add the averages determined 
above to obtain noiseless light curves, with averages that 
satisfy exactly the empirical PLC relation. 

- Add Gaussian noise to the above noiseless light curves. The 
noise realizations used for the light curves with V averages 
are different from the ones used for the light curves with the 
B averages (however, they have the same standard deviation 

Cr). 

We note that the above generation of test light curves is not en- 
tirely consistent, because the noiseless B light curves still have 
Fourier decompositions corresponding to the V light curves. 
However, this inconsistency has only a small effect on the es- 
timated average magnitudes as we have shown in Sect. 5 on a 
smaller set of real B light curves. 

In Fig. 10 we show an example of the improvement ob- 
tained by the application of TFF. The standard deviations 
around the best fitting straight lines are 0.064 and 0.038 for 
the DFF and TFF results, respectively. The derived slopes with 
their l<x standard deviations of the means are the following: 
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-2.485 ± 0.022 (DFF) and -2.490 ± 0.013 (TFF). These slopes 
are within the error limit of the empirical value of -2.467. 
Other realizations with higher noise also show the advantage 
of using TFF over DFF. 

7. Tests on the OGLE LMC data 

To test the applicability of TFF on real astronomical time se- 
ries, we choose the RR Lyrae database of OGLE on LMC 
(Soszynski et al. 2003). Because the database contains more 
than 7000 entries and it is out of the scope of this paper to per- 
form tests on all these stars, we choose a subset of it. This sub- 
set comprises fields # 4, 5 and 6 in the central, well-populated 
region of the LMC bar and fields # 17, 18, 19 and 20 in the far 
less-populated outermost part. Since the observational strategy 
of OGLE preferred the I (Cousins) band, the number of data 
points are different in the various wavebands. The overall num- 
ber of data points, TFF fitting accuracy, QF and the total num- 
ber of objects are (20, 0.07, 30, 1765), (40, 0.06, 35, 1886) and 
(500, 0.06, 78, 1947) in B, V and I colors, respectively. We 
note that the fitting accuracy and QF have substantial star-to- 
star scatter. 
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Fig. 12. Con-elation between Wo = -1.59 log P and W = 
V - 3.1{B - V) (adjusted to field averages) for 218 RRab stars 
of selected LMC OGLE fields. The standard deviation of the 
fit is 0.096 mag, the l<x error of the slope is 0.08. The 45° line 
is shown for reference. The result have been obtained by DFF 
with the top 100 most deviating stars clipped. Filled circles de- 
note variables in field #20. 
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Fig. 11. PLC relations computed from the DFF and TFF anal- 
yses of a subset of the OGLE database of RRab stars in the 
LMC. The selective absorption coefficient Ry is set equal to 
3.1. 



Two tests are performed. First we check how well the PLC 
relations of KW01 can be recovered, then we investigate the 
variation of the Fourier phase of the V light curves with the 
period. Because some light curves suffer from excessive noise 
(or some other types of deformation) we apply a parameter 
filter to keep only the good/reasonable quality light curves. 
In sorting out objects for the check of the PLC relation in 
(B,V) we require: QF > 20.0, cr fit / < 0.02, A x > 0.18, 
0.3 < (B - V)o < 0.4. Here, except for the color index, all 
criteria refer only to the V light curves. The dereddened color 
difference, (B - V)o, is computed from the observed one with 
the assumption of £b-v =0.1 (see, e.g. Kovacs 2000). There 
remained 318 and 334 stars passing these criteria in the DFF 
and TFF analyses, respectively. 



I 




.2 .3 .4 .5 .6 

W + const. 

Fig. 13. Correlation between Wo = -2.33 log P and W = 
V - 3.1(B - V) (adjusted to field averages) for 234 RRab stars 
of selected LMC OGLE fields. The standard deviation of the 
fit is 0.092 mag, the l<x error of the slope is 0.05. The 45° line 
is shown for reference. The result have been obtained by TFF 
with the top 100 most deviating stars clipped. Filled circles de- 
note variables in field #20. 



Figure 1 1 shows the resulting log P — > Wb-v relations for 
these stars. We see that there is no particular improvement in 
the tightness of the relation whether using TFF or DFF. This is 
partially understandable, because, as Table 2 shows, we expect 
only moderate improvements in both colors. Although this im- 
provement is too small to be easily visible in the above plot, 
we can test the difference by employing a direct search for 
the best single-parameter regression. To take into account the 
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possible zero point differences between the different fields, we 
consider each field as a different 'cluster' and employ optimum 
(field-dependent) zero point shifts together with the uniform 
(i.e., field-independent) linear regression. We note that the zero 
point shifts were in general in the range of ±0.05 mag. 
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Fig. 14. Low-order Fourier phases computed by DFF from the 
V light curves of 975 RRab stars of selected LMC OGLE fields. 
All stars have QF > 20 and Ai < 0.18. 
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Fig. 15. Low-order Fourier phases computed by TFF from the 
V fight curves of 1056 RRab stars of selected LMC OGLE 
fields. All stars have QF > 20 and A i < 0.18. 



Figures 12 and 13 show the resulting regressions. We used 
an iterative procedure in which at each step of the iteration we 
discarded the most deviating star from the best regression. We 
repeated this procedure 100 times. Of course, this excessive 
outlier selection is not justified statistically. However, our goal 
is to see which method is capable of estimating the functional 
dependence more accurately. 

We found that in both cases the data quality was good 
enough for the code to select P as the best regression parameter. 
This parameter was always the best among the other Fourier pa- 
rameters at each stage of the iteration. However, as it is shown 
in the figure captions, TFF converged to a more accurate value 
of the slope of the log P — > Wb-v relation. The slope obtained 
by DFF is completely out of range of the expected value of 
-2.47 (see KW01). It is also clear that the formal error substan- 
tially underestimates the true error in this case. The next best 
fitting Fourier parameter yielded standard deviations of 0.11 
and 0.13 mag for the DFF and TFF methods, respectively. If 
we compare these values with the dispersions obtained for the 
regressions with P, we get that the increase is a factor of 1.16 
for DFF, whereas it is 1.37 for TFF. This yields a higher signif- 
icance for the correlation with P. 

It is worthwhile to mention that Soszynski et al. (2003) did 
not derive the log P — > Wb-v relation. However, they did it for 
the (V, I) colors, probably because of the higher accuracy of the 
data in these colors. We note that from the subset used in this 
paper we also derived consistent slopes for the log P — > Wv-i 
relation. We got -2.584 + 0.032 and -2.662 ± 0.031 from the 
DFF and TFF analyses, respectively. The two methods perform 
similarly also in other aspects, but the TFF sample contains 
some 40 more stars, due to the better quality of the TFF fits. 

In both pairs of colors, the resulting regressions display 
larger dispersions than expected from standard statistical es- 
timates. Although there might be several sources of the exces- 
sive scatter, crowding should definitely play a role (see Kiss 
& Bedding 2005). Our lower limit set for the A\ amplitude 
is aimed at filtering out some of the blends in a crude way. 
Obviously, a more sophisticated method is needed to be more 
successful in filtering out blended variables. 

A different test can be performed on the Fourier phases. 
This test is less stringent than the one presented above, be- 
cause the quality of the result is judged from the tightness of 
the phase progression with the period, that is a good criterion 
only if there are reasonable pieces of evidence that the metal- 
licity does not have a large scatter. This assumption is proba- 
bly not a bad one for the LMC (e.g., Gratton et al. 2004 from 
spectroscopy of RRab stars; Kovacs 2001, from double-mode 
stars). In Figs. 14 and 15 we show the progressions obtained 
for the V light curves of the selected fields mentioned at the 
beginning of this section. We see that the TFF data indeed ex- 
hibit a tighter correlation. Furthermore, the original number 
of stars of 1886 reduces in a less extent for TFF when var- 
ious QF cutoff values are used. We also observe some hori- 
zontal structures in the TFF plots (see the sequences of circles 
at <f2\ = 1.95 and 2.15). When lower QF cutoff is used, then 
these (and also some additional) structures become more visi- 
ble. This indicates that the effect is partially due to poor data 
quality. However, visual inspection of the light curves produc- 
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ing these structures, and a similar test performed on the more 
substantial data in I color have shown that there are groups of 
stars producing nearly constant TFF phases even when large 
QF cutoff is used. Furthermore, because of the basic template 
set contains only 248 stars, it is possible that certain special 
structures observed by OGLE can be reproduced only by a few 
stars with very similar Fourier phases (which may yield con- 
stant phases if the template polynomial degree M is equal to 
or 1 - a common case for poor data quality). 

8. Conclusion 

We devised a full-fetched Template Fourier Fitting (TFF) 
method for the computation of the Fourier decompositions of 
undersampled, noisy fundamental mode RR Lyrae (RRab) light 
curves. The method can be extended to other types of vari- 
ables, assuming that the corresponding template set is reason- 
ably complete (i.e. it contains most of the 'flavors' of the po- 
tential targets). The main features of the method are as follows. 

- Finds matches between the target and individual template 
members (i.e., it does not employ multi-template regres- 
sions); 

- Fits templates to the target by applying polynomial trans- 
formation of the template; 

- Optimizes the degree of the polynomial transformation be- 
tween zero and two, depending on the data quality. 

We performed a number of tests to investigate the range of ap- 
plicability of the method. These tests included the estimation 
of: (i) the photometric iron abundance from the ipn Fourier 
phase and period P; (ii) the average magnitudes in various color 
bands; (iii) the period-luminosity-color relation. In all these 
tests TFF proved to perform better than an optimized Direct 
Fourier Fitting (DFF) method when the noise level was high 
or the number of data points was small. The two methods yield 
the same solution for light curves of high Signal-to-Noise Ratio 
(SNR). For example, when estimating [Fe/H] from light curves 
with 30 data points we expect a factor of two increase in the 
accuracy of [Fe/H] if the standard deviation cr of the noise is 
lower than 0.03. For higher noise levels, TFF remains to be 
more accurate up to fairly high number of data points (e.g. with 
cr = 0.06 we get better result with TFF even alN ~ 100). 

When applied to a subset of the OGLE RR Lyrae database 
on the LMC (see Soszynski et al. 2003), TFF produces sta- 
tistically more significant period-luminosity-color relation for 
(B, V) colors, although the significance of the relation still re- 
mains marginal on this subset. For (V, I) colors we get the 
same, statistically significant relation from both DFF and TFF. 
Nevertheless, the dispersion of the relations are high in all col- 
ors, suggesting the importance of crowding effects. 

The stable performance of TFF for undersampled and noisy 
light curves makes it suitable to revisit problems such as the 
RR Lyrae metallicities in globular clusters and in galaxies (e.g. 
in the Magellanic Clouds) or the determination of average col- 
ors and empirical relations. With the various large-scale sur- 
veys (microlensing, variability, transit, etc.) there is an increase 
in the number of the good- and bad-quality light curves alike. 



Therefore, we expect TFF a useful supplementary method to 
the traditional Fourier fitting. 

Acknowledgements. The support of OTKA grant K-60750 is ac- 
knowledged. 
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